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Density distributions of outflow driven turbulence 
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ABSTRACT 

Protostellar jets and outflows are signatures of star formation and promising mech- 
anisms for driving supersonic turbulence in molecular clouds. We quantify outflow- 
driven turbulence through three-dimensional numerical simulations using an isother- 
mal version of the robust total variation diminishing code. We drive turbulence in 
real-space using a simplified spherical outflow model, analyse the data through den- 
sity probability distribution functions (PDF), and investigate the Core Formation 
Rate per free-fall time (CFR//). The real-space turbulence driving method produces 
a negatively skewed density PDF possessing an enhanced tail on the low-density side. 
It deviates from the log-normal distributions typically obtained from Fourier-space 
turbulence driving at low densities, but can provide a good fit at high-densities, par- 
ticularly in terms of mass weighted rather than volume weighted density PDF. Due 
to this fact, we suggest that the CFR// determined from a Fourier-driven turbulence 
model could be comparable to that of our particular real-space driving model, which 
has a ratio of solenoidal to comprcssional components from the resulting turbulence 
velocity fields of ^0.6. 
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1 INTRODUCTION 

Protostellar jets and outflows are important by-products of 
the early stages of star formation where a collapsing proto- 
stellar system expels excess accretion material in order to 
, lose angular momentum. This material takes the form of 
bipolar protostellar jets, which are collimated and acceler- 
ated along the rotation axis of the system. The jets entrain 
ambient cloud material to form protostellar outflows which 
propagate and deposit energy and mom entum up to sev- 
eral p arsecs into the molecular cloud (e.g. iMcGroartv et al.l 
120041 ) . Protostellar jets and outflows are thus promising can- 
didates for turbulence generation in molecular clouds as they 
drive it internally and may enable star formation to be a 
self-regulating process. 

Although simulations of protostellar outflows clearly 
show a transfer o f momentum from the jet to the ambi- 
ent medium (e.g. iMoraghan et al.l 120081 ). previous studies 
have suggested they may n ot be capable of sustai ning turbu- 
lence in molecular clouds. iBaneriee et al.l (|2007i ) simulated 
single jets and analysed the data using velocity probability 
distribution functions (PDF) to find that supersonic fluc- 
tuations decay quickly and do not spread far from the jet. 
They found subsonic non-compressional modes occupy the 
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disturbed volume instead. Later. lCunningham et al.l Q2009J) 
showed that by taking a more realistic non-uniform ambient 
environment consisting of pre-existing turbulent motions or 
'fossil cavities', then the outflow morphology can be signifi- 
cantly modified leading to a more efficient transfer of energy 
to the ambient medium. Recently. ICarroll et alj (|2009, 2010) 
built upon the previous work with simulations of multiple 
interacting outflows in a global parsec scale volume analysed 
via power-spectra. Their simulations suggest that outflows 
can sustain turbulence through their interactions. 

On the observat i onal side, in a recent study by 
lloannidis fc Froebrichl (|2012l ) of the Serpens and Aquila re- 
gions, the authours estimated the outflows they observed 
would be unable to support the turbulence of the surround- 
ing molecular cl ouds. In co n trast, in a study of the Perseus 
cloud complex, lArce et al.l 1)20101 ) found that a significant 
fraction of the turbulence can be sustained by outflows. 
Similarly, Naka mura et al.1 (|2012j ) who surveyed the L1641- 
N cluster suggest that outflow feedback can influence the 
dynamical evolution of the clump. 

Since the work of lVazquez-Semadenil II 19941 ) it is widely 
accepted that for isothermal supersonic turbulence, the den- 
sity PDF is log-normal with a standard deviation, a, pro- 
portional to the Mach number M. The relationship is now 
usually defined throug h the relation a 2 = ln(l + b 2 A4 2 ), 
as first formulated by IPadoan et al.1 |l997) with b — 0.5. 
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Othe r groups find similar re sults but with different b values 
(See iFederrath et al.l IJ2008T ) and references therein). Such 
log-normal fits readily appear, not only in many numeri- 
cal simulations of turbulence, but also through analytical 
analysis, and observational data of molecular clouds where 
observers measure column -densities of tracer molecules (e.g. 
iFroebrich fc Rowlesllioiol) . 

Many recent works on numerical simulations of tur- 
bulence now drive it in Fourier-space through a combina- 
tion of a s olenoidal and compressional based forcing mech- 



anism c 



ilcnoidal and compressional based forcing mech - 
Federrath et~"al"1l200Sl . I2OI0I ; iMolina et al.ll2012l ). 



IFederrath et al.1 (J2008 ) find that the addition of compres- 
sive forcing makes the standard-deviation of the density 
PDF about 3 times wider than that of solenoidal driving 
alone, and also confirm the relationship between standard- 
deviation and Mach number. However, little research has 
been performed relating to the density PDF from simula- 
tions of turbulence driven in real-space. 

In this work we investigate the concept of outflow driven 
turbulence numerically using a modified version of the 
isothermal total variation diminishing (tvd) code. Whereas 
most previous turbulence simulations drive the turbulence 
in Fourier-space, we use a simple model considering proto- 
stellar outflows to drive the turbulence real-space. We char- 
acterise the resulting turbulence in terms of density PDF 
and core formation rate per free fall time (CFR//). We find 
that driving the turbulence in real-space via outflows, leads 
to density PDFs with negatively skewed log-normal fits pos- 
sessing enhanced tails at the low-density side. This devi- 
ates from the log-normal distributions commonly found in 
isothermal driven turbulence simulations. However, due to 
good agreement on the high-density side of the density PDF, 
particularly in terms of mass, the results suggest that esti- 
mates of the CFR// based purely on log-normal fits, or from 
Fourier driven turbulence, could still provide comparable re- 
sults to those of real-space driven turbulence. 

The layout of our paper is as follows. We describe the 
basis of our model for generating turbulence and the com- 
putational code used in this study in § [2] The numerical 
results and density PDFs are presented and discussed in 
§ [3] We summarise and conclude our findings in § [4] 



2 NUMERICAL METHODS 

Observations of many active molecular clouds, such as the 
Serpens and Aquila region, have revealed many complex dy- 
namical processes including starless cores, young stellar ob- 
jects , outflows, and turbulence (e.g. Iloannidis fc Froebrichl 
|2012j). Modelling a realistic molecular cloud would be very 
comprehensive, requiring an external medium to sufficiently 
bind it, as well as magnetic fields, molecular physics, self- 
gravity, and a star formation accretion-ejection mechanism. 
As a first step, we study a simplified case as we are 
trying to understand the properties of turbulence driven by, 
and associated with, outflows only. We therefore inject a se- 
ries of spherical outflows without self-gravity or magnetic 
field effects into the computatio nal domain. As such, sp her- 



ical outflows w ere consider ed by Li fc Nakamural (|2006l ) and 
iMatznerl 1120071 ). Although IMatznerl (120071 ) claim collimated 
outflows are more efficient than spherical outflows at driving 
turbulence as they can propagate further from their source 



thus driving turbulence on larger scales which decay slower, 
in this work our purpose is to start with a simple model. 

Our turbulence forcing model proceeds as follows; A 3D 
computational domain is initialised with a uniform gaseous 
medium. Different points within the domain are randomly 
selected. At each randomly selected point, the mean density, 
p, of a volume defined by the outflow radius, R, is deter- 
mined. If the mean density of this volume is larger than the 
global mean, po, and no point within the volume is below 
a lower-threshold value O.OOlpo, then it is chosen to be the 
location of an outflow. 

The momentum, P, which an outflow should introduce 
is defined as follows, 



P(r)vj(R) (^Y dV 



(1) 



where p(r) is the density within the outflow volume depend- 
ing on distance vector r from the outflow center, and Vj (R) 
is the outflow velocity at its outermost boundary, R. The 
q parameter, set to 1.0, attempts to mimic the observed 
'Hub ble-law' type vel ocity observed in protostellar outflows 
(e.g. lArce et alj|2007l ). 

As the total scalar momentum, P, which each outflow 
should introduce is known, vj (R) can be expressed in terms 
of P and the integral part of Eq.fT] However, by setting a 
new velocity field, the total momentum of the outflow may 
not be conserved due to a non-uniform density distribution 
existing within the chosen outflow radius, especially dur- 
ing later times of the simulation. Therefore we calculate the 
added momentum and then subtract it within the radius R. 
This ensures the vector sum of the momentum is zero and 
the total added momentum on the computational domain 
remains conserved. 

In order to give physical meaning to our model, we make 
use of the s ame dimensionle ss parameter sche me as first in- 
trodu c ed by IMatznerl (|2007l ) and later used bv lCarroll et al.l 
(2009. 120101 ). Following this scheme conveniently links our 
results to physical quantities and enables comparison of our 
results to those of different codes. We thus define in our 
model the dimensionless quantities of mass, mo, length, Zo, 
and time, to, as, 



m 
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/<) 



p- 
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and to 



P*S? 



(2) 



We assign the values of density p = 2.51 x 10 -20 gem -3 , 
outflow momentum P — 3.98 x 10 39 gem -3 s~ , and outflow 
rate per u nit volume S = 6.31 x 10~ cm - s _1 , which were 
chosen bv lCarroll et al.l (|2010l ) to approximate the physical 
values in a typical star forming region. After implementing 
the above unit scheme, the resulting physical unit quantities 
that we use in our code are; mo = 18.7 Mq, lo — 0.37 pc, to 
- 0.34 Myr, and Vo - lo/to = 1.064 kms -1 . The local sound 
speed has a value of 0.587 kms -1 . We set our domain scale 
length as 8I0 (3.96 pc), choose an outflow radius to be R = 
0.3 Zo, and run the simulation for time 3io (~1.02 Myr). 

We perform the simulations on a fixed uniform grid of 
512 3 cells with periodic boundary conditions to emulate a 
subset of a molecular cloud. Numerical convergence of the 
density PDF data was tested at lower resolutions of 128 3 
and 256 3 . Although the majority of the cells provide precise 
convergence forming the peaks of the density PDF, there 
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is a slight dis crepancy in the t ails. T he same phenomenon 
was noted by iFederrath et al.l (120101 ) who decides it is dif- 
ficult to obtain exact numerical convergence through such 
turbulence simulations. However, we conclude a resolution 
of 512 3 is sufficient for our purposes as both the fitted mean 
and standard deviations of the PDF at all resolutions agree 
to four decimal places. 

Finally, the model was implemented using an isother- 
mal version of th e tvd multi-dimensional fixed-grid code by 
iKim et al.l |l999|). The code is based on an explicit finite- 
difference scheme which uses a second-order accurate Roe- 
type up-winded Riemann solver to calculate the inter-cell 
flux. The performance of the code has been compared with 
several oth er prominent hydrody namic codes in use today in 
a study bv lKitsionas et al.l ((2009J). The study included three 
of the major grid-based codes (enzo, flash and zeus) and 
compared simulations of decaying, isothermal, supersonic 
turbulence using the same initial conditions. In all cases, 
the tvd code was found to give comparable results to its 
contemporaries but with substantially shorter run-times. 
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Figure 1. Variation of the root mean square Mach number, 
M, with time during our simulation. The outflows quickly gener- 
ate turbulent motions which saturate after a time of 0.5 to- The 
range over which we produce our time-averaged density PDFs is 
highlighted between the vertical dashed lines between time 1.0 
to 3.0to- The horizontal dotted line represents the Mach number 
corresponding to the outflow velocity scale, vq. 



3 RESULTS 



The evolution of the root mean square Mach number, Ml, 
during the 3.0 to running time is shown in Fig.Q] The high 
velocity spikes visible in Fig.Q] are a consequence of our 
method of implementing the outflows. An outflow radius 
may include some cells of very low-density, and through the 
conservation of momentum these cells would have very high- 
velocity. We can control the magnitude of the high- velocity 
spikes by increasing the lower-threshold value where an out- 
flow can be placed, but we compromise to ensure the out- 
flows may be placed as randomly as possible. The very high- 
velocity cells are localised and have little effect on the global 
turbulence evolution. The outflows themselves are also quite 
localised so although normally possessing an initially high 
Mach number of M ~ 20, the global average M is much 
lower, yet still supersonic at ~1.7. 

A volume rendering of density during our simulation 
is shown at a time of 2.0 to in Fig. [2] It is apparent that 
the majority of the volume consists of gas at low-density 
with some spherical outflow cavities visible. To analyse 
the data quantitatively, we plot normalised Volume and 
Mass weighted density PDFs in natural-log. Fig.[3K is the 
time-averaged Volume weighted density PDF (fronr 25 data 
dunrps) between a simulation time of 1.0-3.0 to- Fig.^ is 
the corresponding Mass weighted density PDF. The solid 
black lines represent the density PDF of the normalised 
data arrays, p s (s), where s — ln(p/po). We measure the 
Mean, p, — s — J sp s (s) ds, and Standard Deviation, a 2 
— J (s — s) p s (s) ds of the distributions. We also deter- 
mine the best log-normal fit to the distributions using the 
Levenberg-Marquardt method. The resulting log-normal fits 
are plotted as the red curves. A log-normal distribution be- 
ing defined as: 



f(s) 



v^- 



zexp 



— (s — pln) 
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(3) 



where p^x and a\ N are the Mean and Standard Deviations 
of the log-normal fit. 



ft' 1 * 1N "L 



mi 



...?■»' 







Figure 2. A 3D volume rendering of our simulation depict- 
ing density at time 2.0 to. Despite a large density contrast from 
2.55x 10~ 5 to 61.72 po, the majority of the volume is occupied by 
low-density spherical cavities interspersed by filamentary struc- 
tures. 



Also determined are the values of Skewness and Kur- 
tosis. Skewness is a measure of symmetry of the distri- 
bution and defined as <S S = J [(s — s) /a] p s (s) ds. Kur- 
tosis measures the extent to which a distribution devi- 
ates from a normal-distribution and is defined as /C s = 
J [(s — s) /a] p a (s) ds — 3, where the last term is often in- 
cluded to ensure the Kurtosis of a standard normal distri- 
bution is 0. The measured quantities are listed in Table [T] 

The negative skewness values indicate the density PDFs 
are skewed to the low-density side, although the skewness of 
the Mass weighted density PDF is closer to log-normal. This 
could be expected due to the mass of the diffuse material 
having less effect on the low-density side of the density PDF. 
The Kurtosis measurements show the same trend while be- 
ing positive or 'leptokurtic' which indicate an enhancement 
in the density PDF tails. 
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Figure 3. (a) Normalised density PDF plot in log. Black line 
represents the time-averaged Volume fraction. Red curve repre- 
sents the best log-normal fit to the data using the Levenberg- 
Marquardt method. Blue curve is a log-normal approximation 
assuming Fourier-based driving at the same Mach number. Dot- 
ted dark-blue and dashed light-blue curves are the log-normal 
approximations of purely compressible and solenoidal driven tur- 
bulence respectively, (b) Same as (a) but for Mass fraction with 
a vertical green line to represent the threshold density for the 
CFR// measurement. 



Comparing the density PDFs to their log-normal fits in 
Fig. [3] we see a distinctive low-density tail on the low-density 
or diffuse side, pronounced by the log-based plot. It is due 
to the spherical outflows of our model which have expanded 
outwards leaving lower density 'fossil cavities'. In our simu- 
lations we see how fossil cavities do enhance the transfer of 
momentum to the ambient medium. Through the conserva- 
tion of momentum, material moves faster in the low-density 
regions, thus enabling the outflow shocks to propagate fur- 
ther from the source through older cavities than they would 
otherwise travel in an undisturbed uniform medium. 

Due to the fact we do not implement self-gravity, 
we do not obtain an enhanced high-d ensity tail on the 
PDF as found by o t her a uthours (e.g. ICho fc Kirnl l201ll ; 
iFederrath fc Klessenl l2015t ). Although we find the Mass 
weighted density PDF is a better match to the high-density 
side of the log-normal fit than the Volume weighted equiva- 
lent. Again, this is due to the Volume weighted density PDF 
being skewed towards lower density by a greater amount. 

We also plot a log-normal curve to represent the re- 
sult one would obtain from Fourier driven turbulence at the 
same average Mach number as our simulation. We proceed 
by makin g use of the forcing pa rameter, £, originally em- 
ployed by IFederrath et al.l (|2008l ) as the relative strength of 
the solenoidal to compressional components. Although the 
authors defined it with respect to their forcing function, we 
define it at the end of our simulation based on the result- 
ing turbulent velocity field averaged over the time interval 
1.0-3.0 to. We measure £ to be 0.614, a value that indicates 



solenoidal modes slightly dominate in our fully converged 
sa turated turbulent velo city field. Then through equation 5 
of lFederrath et al.l ((2008) we determine the corresponding b 
value (0.591) and use it, along with a measure of our average 
Mach number (1.69), in the a 2 = In (l + b 2 M 2 ) equation, 
and th e \i — —a 2 /2 relation of IPassot fc Vazquez-Semadenil 
(1998), to obtain a standard-deviation and mean value for 
plotting a general Fourier- driven based log-normal approxi- 
mation. In Fig.|3K we also plot the log-normal curves cor- 
responding to the two extremes; pure solenoidal driving 
(£=1.0) and pure compressible driving (£=0.0). 

From visual inspection, we see the Fourier-driven based 
log-normal approximation is a more precise fit to the high- 
density end of the Volume weighted density PDF than its 
best log-normal fit. In the Mass weighted density PDF, the 
Fourier-driven based log-normal approximation is a closer 
match to the best log-normal fit (see standard deviation val- 
ues in Table[T|. 

The departure of our data from the standard log-normal 
shaped density PDF is most likely due to the fact that we 
drive the turbulence in real-space using outflows, whereas 
most numerical works who study density PDFs drive the 
turbulenc e in Fourier-space based on Gaussian distribu- 
tions (e.g. IFederrath et alll2008l ; IFederrath fc Klessenl I20I2I ; 



iMolina et al.ll2012f ). By driving turbulence in Fourier-space, 
the material over the entire domain is mixed in a sinusoidal 
based manner; some portion becomes high-density, and an 
equal portion becomes low-density. It is hard to obtain a 
significantly non log-normal density PDF in this way with- 
out extra physics such as self-gravity. We would expect a 
difference by driving the turbulence in real-space, as our 
spherical outflow model is highly localised. The outflows are 
introduced to regions of above average density and easily 
inflate low-density bubbles in the computational domain. 

The relationship between the star formation rate and 
density PDF has been inve s tigate d by several groups, 



namel y, IKrumholz fc McKeel (120051). iPadoan fc Nordlundl 
j|201ll ) and lHennebelle fc Chabried (201ll ). Each group used 
a slightly different theoretical measure to define the criti- 
cal density for star formation to occur and thus estimates 
different values for the star formation rate. 

Here we investigate the theory of CFR//. In our Mass 
weighted density PDF plot of Fig.[3j3 we choose a criti- 
cal density of p/po = 10 and calculate the fraction of the 
PDF area greater than this density. This is equivalent of 
imple menting the Error Functio n as defined by equation 
20 of IKrumholz fc McKeel (|2005T ) with e core = 1, and later 
used bv lCho fc Kim (2011). Although our model is not en- 
tirely appropriate for a precise estimate of the CFR// as 
we do not implement self-gravity and so can not take a 
critical density from the Je ans Length measurements as in 
IKrumholz fc McKeel (J2005I ). here we simply take a critical 
density independent of Mach number and can still make 
a comparison of this arbitrary CFR// between the Mass 
weighted density PDF and its corresponding log-normal fit- 
tings. 

The CFR// values are listed in Tabled] The results 
show that CFR// of the Fourier driven log-normal fit based 
on the £ value determined from our resulting velocity fields, 
differs from that of the outflow driven data by a factor of 
~2. But as we are using a slightly different definition of 
£ as explained above, a factor of 2 difference may still be 
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Table 1. Summary of measurements obtained from the den- 
sity PDFs of the outflow driven (OD) data, the best log-normal 
(LN) fits, and the Fourier-driven based log-normal approximation 
(FD), for both the Volume and Mass weighted cases. Bottom row 
shows the estimated values of the CFR// , which is a measurement 
based on mass. 





(a) Volume weighted 


(b): 


Mass weij 


jhted 




OD LN 


FD 


OD 


LN 


FD 


Mean /i 


-0.732 -0.272 


-0.348 


0.501 


0.539 


0.348 


Std-dev a 


1.560 1.057 


0.834 


0.896 


0.845 


0.834 


Skewness S a 


-1.490 0.0 


0.0 


-0.326 


0.0 


0.0 


Kurtosis Kk 


3.016 0.0 


0.0 


0.99 


0.0 


0.0 


CFR// 






0.0185 


0.0189 


0.0097 



in good agreement. Overall, a log-normal approximation to 
the data, or a Fourier-driven turbulence model could still 
provide comparable CFR// results with our particular real- 
space driving model as the CFR// is dependent on mass 
measurements at high-densities. We will be able to perform 
a more accurate analysis in a future publication with self- 
gravity implemented. 



include the effects of magnetic fields as are present in molec- 
ular clouds. 
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4 CONCLUSIONS 

We present simulations of outflow driven turbulence per- 
formed using the isothermal tvd code. The turbulence is 
driven exclusively in real-space by a series of spheri cal out- 
flows based on the parameter scheme as used by iMatznerl 
i|2007h and ICarroll et all (|2009l , I2010T ) in order to approxi- 
mate the physical values of a typical of star-formation cloud. 

We analyse the resulting turbulence in terms of Volume 
and Mass weighted density PDFs, and the CFR// over a 
driven regime in which outflows are continuously introduced. 
Our primary finding is a departure from log-normality of the 
density PDF due to an extended low-density tail, quantified 
by measures of negative Skewness and positive Kurtosis. It 
is a consequence of our localised real-space driving method 
where the multiple outflows sweep the ambient gas into thin 
shells, thus enhancing the volume of low-density cavities. 
We propose that if there is localised turbulent driving in 
an active cloud via an outflow, stellar wind, or supernova, 
then the density PDF may have an enhanced tail at the 
low-density side. However, log-normal fits to the data may 
still provide good approximations to the high-density side 
especially when considering Mass weighted density PDFs 
rather than Volume weighted. 

Although we do not implement self-gravity, we quan- 
tify our results with a CFR// calculation which we define 
as the fraction of density PDF area greater than p /pp = 10, 
as similarly implemented by ICho fc Kiml (|201ll ). Our re- 
sults suggest that the CFR// based on our particular simple 
real-space driving model with £ ~ 0.6 is comparable to the 
CFR// that would be obtained through Fourier-driving at 
the same Mach number. 

In future work, we shall perform a power-spectra anal- 
ysis of higher-resolution simulations and a more detailed 
analysis of the compressional and solenoidal generated by 
our real-space driving mechanism. We will also improve our 
physical model by implementing self-gravity so that dense 
cores form naturally, add collimation to the outflows, and 
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